---
title: 'Replication script for: Taxing Your Cake and Growing It Too:  Public Beliefs
  on the Dual Benefits of Progressive Taxation'
author: "Bastian Becker, Bruno Castanho Silva, Hanna Lierse"
date: "`r Sys.Date()`"
output: 
 html_document:
   theme: flatly
   highlight: tango
   toc: true
   toc_float: true
   code_folding: show
   number_sections: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = F, warning = F)
```

# Packages and data

If you have this file saved in the same folder as the dataset file `replication_data_jpp.csv`, all commands below will run and produce the respective results.

The packages necessary for replicating this analysis are:

```{r}
library(tidyverse)
library(cregg)
library(dotwhisker)
library(gridExtra)
```

Data:

```{r}
data_conjoint <- read.csv('replication_data_jpp.csv', stringsAsFactors = T)
```

## Data transformation

Create one dataset where the reference category is tax increases, and one where it is decreases, to make the estimation and plotting easier:

```{r}
data_conjoint_increases<-data_conjoint
data_conjoint_increases$einkom<-relevel(data_conjoint_increases$einkom,"Income - General Increase")
data_conjoint_increases$untern<-relevel(data_conjoint_increases$untern,"Corporate - General Increase")
data_conjoint_increases$erbsch<-relevel(data_conjoint_increases$erbsch,"Inheritance - General Increase")

data_conjoint_decreases<-data_conjoint
data_conjoint_decreases$einkom<-relevel(data_conjoint_decreases$einkom,"Income - General Decrease")
data_conjoint_decreases$untern<-relevel(data_conjoint_decreases$untern,"Corporate - General Decrease")
data_conjoint_decreases$erbsch<-relevel(data_conjoint_decreases$erbsch,"Inheritance - General Decrease")
```

# Figure 1

Estimate results for both growth and equality:
```{r}
# Estimate main results
# Tax increases
results_main <- cj(reduneq ~ 
                     einkom +
                     erbsch + 
                     untern, 
                   data = data_conjoint_increases, id = ~id, estimate = "amce")
names(results_main)[4]="term"

results_growth <- cj(incgrow ~ 
                       einkom +
                       erbsch + 
                       untern, 
                     data = data_conjoint_increases, id = ~id, estimate = "amce")
names(results_growth)[4]="term"

results_increases <- bind_rows(results_main, results_growth) %>% subset(grepl("Top Increase", .[,"term"]))
# Tax decreases
results_main <- cj(reduneq ~ 
                     einkom +
                     erbsch + 
                     untern, 
                   data = data_conjoint_decreases, id = ~id, estimate = "amce")
names(results_main)[4]="term"

results_growth <- cj(incgrow ~ 
                       einkom +
                       erbsch + 
                       untern, 
                     data = data_conjoint_decreases, id = ~id, estimate = "amce")
names(results_growth)[4]="term"

results_decreases <- bind_rows(results_main, results_growth) %>% subset(grepl("Top Decrease", .[,"term"]))
# Plot results
results_both <- bind_rows(results_increases, results_decreases)

brackets <- list(c("Inheritance Tax", "Inheritance - Top Increase", "Inheritance - Top Decrease"),
                 c("Income Tax", "Income - Top Increase", "Income - Top Decrease"),
                 c("Corporate Tax", "Corporate - Top Increase", "Corporate - Top Decrease"))

names(results_both)[1]='model'

## Invert levels of increase/decrease:
results_both$term <- as.factor(as.character(results_both$term))
results_both$term <- ordered(results_both$term, levels = c("Income - Top Increase","Income - Top Decrease","Inheritance - Top Increase","Inheritance - Top Decrease",
                                                           "Corporate - Top Increase","Corporate - Top Decrease"))
```

Plot:

```{r}
p.main <- dwplot(results_both, ci=.84, dot_args = list(size = 2, aes(color=model))) + theme_bw() +
  geom_vline(xintercept=0, colour = "gray20", linetype = 'dashed') +
  xlab("In relation to general tax increase/decrease, \nprobability that will promote...") + ylab("") +
  ggtitle("") + scale_color_manual(values = c("black",'gray50'),
                                   labels = c('...growth','...equality')) +
  theme_minimal() + 
  theme(plot.title = element_text(face="bold"),
        legend.position = 'bottom', legend.title = element_blank(),
        axis.text=element_text(colour="black")) + 
  scale_y_discrete(labels = c('Income - Top Increase'='Top Increase','Inheritance - Top Increase'='Top Increase',
                              'Corporate - Top Increase'='Top Increase','Income - Top Decrease'='Top Decrease',
                              'Inheritance - Top Decrease'='Top Decrease','Corporate - Top Decrease'='Top Decrease')) 


p.main <- add_brackets(p.main, brackets)

p.main
```


# Figure 2

```{r}
# Estimate subgroup results, by ideology
# Tax increases
results_main <- cj(reduneq ~ 
                     einkom +
                     erbsch + 
                     untern, 
                   data = data_conjoint_increases, id = ~id, estimate = "amce", by=~lr.r)
names(results_main)[5]="term"

results_growth <- cj(incgrow ~ 
                       einkom +
                       erbsch + 
                       untern, 
                     data = data_conjoint_increases, id = ~id, estimate = "amce", by=~lr.r)
names(results_growth)[5]="term"

results_increases <- bind_rows(results_main, results_growth) %>% subset(grepl("Top Increase", .[,"term"]))
# Tax decreases
results_main <- cj(reduneq ~ 
                     einkom +
                     erbsch + 
                     untern, 
                   data = data_conjoint_decreases, id = ~id, estimate = "amce", by=~lr.r)
names(results_main)[5]="term"

results_growth <- cj(incgrow ~ 
                       einkom +
                       erbsch + 
                       untern, 
                     data = data_conjoint_decreases, id = ~id, estimate = "amce", by=~lr.r)
names(results_growth)[5]="term"

results_decreases <- bind_rows(results_main, results_growth) %>% subset(grepl("Top Decrease", .[,"term"]))
# Plot results
results_both <- bind_rows(results_increases, results_decreases)


names(results_both)[2]='model'

results_both$term <- as.factor(as.character(results_both$term))
results_both$term <- ordered(results_both$term, levels = c("Income - Top Increase","Income - Top Decrease",
                                                           "Inheritance - Top Increase","Inheritance - Top Decrease",
                                                           "Corporate - Top Increase","Corporate - Top Decrease"
))

results_both <- mutate(results_both, lr.r = paste0(lr.r, ' Respondents'))
```

Plot:

```{r}
p.ideo2 <- mutate(results_both, model = case_when(model == 'reduneq' ~ '... equality',
                                                  model == 'incgrow' ~ '... growth')) %>%
  ggplot(aes(x = estimate, y = term, color = BY)) + facet_wrap(. ~ model) + 
  geom_vline(aes(xintercept = 0), linetype = 'dashed') +
  geom_pointrange(aes(xmin = estimate - 1.41*std.error, xmax = estimate + 1.41*std.error),
                  position=position_dodge(.4)) + 
  theme_minimal() + scale_color_manual(values = c("black",'gray50')) +
  ggtitle("In relation to general tax increase/decrease, \nprobability that will promote...") + 
  ylab(NULL) + xlab(NULL) +
  theme(plot.title = element_text(face="bold"),
        legend.position = 'bottom', legend.title = element_blank(),
        axis.text=element_text(colour="black")) + 
scale_y_discrete(limits = rev)
p.ideo2
```

# Figure 3

```{r}
## Income:
# Tax increases
results_main <- cj(reduneq ~ 
                     einkom,
                   data = data_conjoint_increases, id = ~id, estimate = "amce", by=~high_inc)
names(results_main)[5]="term"

results_growth <- cj(incgrow ~ 
                       einkom, 
                     data = data_conjoint_increases, id = ~id, estimate = "amce", by=~high_inc)
names(results_growth)[5]="term"

results_increases.inc <- bind_rows(results_main, results_growth) %>% subset(grepl("Top Increase", .[,"term"]))
# Tax decreases
results_main <- cj(reduneq ~ 
                     einkom , 
                   data = data_conjoint_decreases, id = ~id, estimate = "amce", by=~high_inc)
names(results_main)[5]="term"

results_growth <- cj(incgrow ~ 
                       einkom, 
                     data = data_conjoint_decreases, id = ~id, estimate = "amce", by=~high_inc)
names(results_growth)[5]="term"

results_decreases.inc <- bind_rows(results_main, results_growth) %>% subset(grepl("Top Decrease", .[,"term"]))

# Bind:
results_inc <- bind_rows(results_increases.inc, results_decreases.inc)

names(results_both)[2]='model'

## Inheritance:
# Tax increases
results_main <- cj(reduneq ~ 
                     erbsch,
                   data = data_conjoint_increases, id = ~id, estimate = "amce", by=~high_wealth)
names(results_main)[5]="term"

results_growth <- cj(incgrow ~ 
                       erbsch, 
                     data = data_conjoint_increases, id = ~id, estimate = "amce", by=~high_wealth)
names(results_growth)[5]="term"

results_increases.inh <- bind_rows(results_main, results_growth) %>% subset(grepl("Top Increase", .[,"term"]))
# Tax decreases
results_main <- cj(reduneq ~ 
                     erbsch , 
                   data = data_conjoint_decreases, id = ~id, estimate = "amce", by=~high_wealth)
names(results_main)[5]="term"

results_growth <- cj(incgrow ~ 
                       erbsch, 
                     data = data_conjoint_decreases, id = ~id, estimate = "amce", by=~high_wealth)
names(results_growth)[5]="term"

results_decreases.inh <- bind_rows(results_main, results_growth) %>% subset(grepl("Top Decrease", .[,"term"]))

# Bind:
results_inh <- bind_rows(results_increases.inh, results_decreases.inh)

## Companies:
# Tax increases
results_main <- cj(reduneq ~ 
                     untern,
                   data = data_conjoint_increases, id = ~id, estimate = "amce", by=~large_comp)
names(results_main)[5]="term"

results_growth <- cj(incgrow ~ 
                       untern, 
                     data = data_conjoint_increases, id = ~id, estimate = "amce", by=~large_comp)
names(results_growth)[5]="term"

results_increases.com <- bind_rows(results_main, results_growth) %>% subset(grepl("Top Increase", .[,"term"]))
# Tax decreases
results_main <- cj(reduneq ~ 
                     untern , 
                   data = data_conjoint_decreases, id = ~id, estimate = "amce", by=~large_comp)
names(results_main)[5]="term"

results_growth <- cj(incgrow ~ 
                       untern, 
                     data = data_conjoint_decreases, id = ~id, estimate = "amce", by=~large_comp)
names(results_growth)[5]="term"

results_decreases.com <- bind_rows(results_main, results_growth) %>% subset(grepl("Top Decrease", .[,"term"]))

# Bind:
results_com <- bind_rows(results_increases.com, results_decreases.com)


## Plot 3 panels:
results <- bind_rows(results_inc, results_inh, results_com)

names(results)[2]='model'

results$term <- as.factor(as.character(results$term))
results$term <- ordered(results$term, levels = c("Income - Top Increase","Income - Top Decrease",
                                                           "Inheritance - Top Increase","Inheritance - Top Decrease",
                                                           "Corporate - Top Increase","Corporate - Top Decrease"
))
```

Plot:

```{r}
p.self2 <- mutate(results, model = case_when(model == 'reduneq' ~ '... equality',
                                                  model == 'incgrow' ~ '... growth')) %>%
  ggplot(aes(x = estimate, y = term, color = BY)) + facet_wrap(. ~ model) + 
  geom_vline(aes(xintercept = 0), linetype = 'dashed') +
  geom_pointrange(aes(xmin = estimate - 1.41*std.error, xmax = estimate + 1.41*std.error),
                  position=position_dodge(.4)) + 
  theme_minimal() + scale_color_manual(values = c("black",'gray50')) +
  ggtitle("In relation to general tax increase/decrease, \nprobability that will promote...") + 
  ylab(NULL) + xlab(NULL) +
  theme(plot.title = element_text(face="bold"),
        legend.position = 'bottom', legend.title = element_blank(),
        axis.text=element_text(colour="black")) + 
  scale_y_discrete(limits = rev)
p.self2
```